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ABSTRACT 

In  this  work,  force  fields  are  developed  for  four 
compounds  categorized  by  the  Army  as  potential 
insensitive  munitions  compounds,  which  includes,  2,4- 
dinitroanisole  (DNAN),  N-methyl-p-nitroaniline  (MNA), 
Dinitropyrazole  (DNP)  and  Nitrotriazolone  (NTO). 
These  force  fields  are  used  to  predict  pharmacokinetic 
properties  such  as  octanol-water  partition  coefficients 
and  Henry's  law  constants  along  with  other 
thermophysical  properties  such  as  vapor-liquid 
equilibria,  vapor  pressure,  critical  parameters  and  normal 
boiling  points. 


1.  INTRODUCTION 

Over  the  past  few  years,  there  has  been  considerable 
interest  in  the  development  of  new  Insensitive  Munitions 
(IM)  that  are  highly  efficient  and  relatively  safe  to 
operate  with  low  chances  of  an  accidental  initiation. 
They  possess  higher  thermal  stability  and  lower  shock 
sensitivity  than  traditionally  used  explosive  compounds, 
such  as  TNT.  Although  insensitive  to  shock,  impact  or 
thermal  effects,  they  are  still  known  to  perform  well.  Due 
to  increased  environmental  and  safety  concerns,  it  is 
necessary  to  determine  how  these  compounds  behave  in 
the  environment  in  terms  of  bioaccumulation  potential 
and  aqueous  solubility  before  being  extensively 
employed  in  explosive  formulations. 

Gillett  has  proposed  a  “pre -biological  screen”,  which 
correlates  the  Henry’s  law  constant  (log  H ),  octanol- 
water  partition  coefficient  (log  Kow)  and  the  material 
half-life  to  predict  which  materials  have  the  potential  to 
be  problematic  if  released  to  the  environment  (Gillett, 
1984).  A  wide  variety  of  experimental  techniques 
(Sangster,  1997)  exist  to  measure  the  partition 
coefficients.  An  alternative  strategy  is  discussed  in  this 
work:  the  application  of  atomistic  computer  simulation  to 
the  prediction  of  partition  coefficients  and  phase 
behavior.  Computer  simulation  has  a  few  benefits  over 
experiments.  Given  an  accurate  intermolecular  potential 


function,  it  is  possible  to  predict  the  behavior  of  a  material 
in  the  environment  before  it  has  been  synthesized, 
allowing  for  a  prescreening  of  potential  candidate 
molecules.  Such  a  prescreening  is  expected  to  lead  to  cost 
savings  by  reducing  the  pool  of  candidate  molecules  for 
synthesis,  and  potentially  leading  to  reductions  in  costs  due 
to  remediation  efforts  near  load-and-pack  facilities. 

In  this  era  of  supercomputing,  complex  processes  can 
be  simulated  in  reasonable  time  scales  with  high  precision 
and  accuracy  across  a  parallel  platform  of  multiple 
processors  with  an  accurate  molecular  model  or  force  field 
and  an  appropriate  simulation  methodology. 

The  structures  of  the  compounds  of  interest  are  shown  in 
Fig.  1. 


Fig.  1.  Structure  of  the  energetic  materials,  a)  DNAN;  b)  MNA; 
c)  DNP;  d)  NTO. 


2.  FORCE  FIELD  DEVELOPMENT 

The  accuracy  of  molecular  simulation  is  dependent  entirely 
on  the  force  fields  used  to  describe  the  interactions  of 
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atoms  with  each  other.  These  interactions  can  be  split 
into  two  categories:  bonded  and  non-bonded.  Bonded 
interactions  account  for  the  conformational  structure  of 
the  molecule,  and  include  bond  stretching,  bond  bending 
and  torsional  rotation  around  the  various  bonds.  These 
interactions  are  typically  described  by  harmonic  or 
cosine  series  potentials,  with  parameters  fit  to  reproduce 
quantum  chemical  data,  such  as  vibrational  frequencies 
and  rotational  barriers.  Non-bonded  interactions  include 
van  der  Waals  and  electrostatic,  which  are  represented  in 
this  work  by  Lennard-Jones  potentials  combined  with 
point  charges.  While  bonded  parameters  are  generally  fit 
to  quantum  chemical  data,  non-bonded  interactions  are 
typically  fit  to  reproduce  experimentally  measured 
quantities,  such  as  liquid  densities  and  heats  of 
vaporization. 

2.1  Non-bonded  interactions 


Non-bonded  interactions  between  atoms  in  each 
molecule  were  represented  with  a  standard  12-6 
Lennard-Jones  potential  (Proceedings  of  the  physical 
society,  1931)  with  a  coulombic  term  for  partial  charges. 
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where  r,j,  sy,  a,,,  q,.  q,  are  the  atom-atom  separation,  LJ 
well  depth,  LJ  diameter  and  the  partial  charges  on  the 
atoms  i  and  j,  respectively.  e0  is  the  permittivity  of 
vacuum.  Initial  estimates  of  the  partial  charges  for  each 
molecule  were  determined  through  a  CHELPG  analysis 
(Breneman  et  al.,  1990)  of  ab  initio  calculations 
performed  at  the  HF/6-3  lg+(d,p)  level  of  theory  and 
basis  set.  This  methodology  was  chosen  based  on  prior 
calculations  for  polar  molecules  such  as  H2S  [Kamath  et 
al.,  2005]  and  dimethyl  ether  [Ketko  et  al.,  2007]. 
During  the  parameter  optimization  process,  the  quantum 
mechanically  derived  partial  charges  were  rescaled  by  a 
factor  of  0.94  to  improve  the  reproduction  of 
experimentally  determined  octanol-water  partition 
coefficients.  Hydrogen  atoms  bonded  to  carbon  atoms 
were  grouped  with  carbon  atoms  to  form  a  single 
“united-atom”.  Lennard-Jones  parameters  a  and  e  for 
each  interaction  site  were  transferred  from  analogous 
compounds  parameterized  previously  in  the  development 
of  the  TraPPE  force  field  (Martin  et  al.,  1998;  Wick  et 
al.,  2000;  Stubbs  et  al.,  2004;  Wick  et  al.,  2005). 


Cross  interaction  parameters  for  unlike  atoms  were 
determined  through  Lorentz-Berthelot  combining  rules 
(Lorentz,  1881;  Berthelot,  1898). 
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Equilibrium  bond  lengths  and  bond  angles,  when  available, 
were  taken  from  the  TraPPE  force  field.  Any  missing 
equilibrium  bond  lengths  or  angles  were  determined  from 
ab  initio  HF/6-3  lg+(d,p)  calculations.  A  harmonic  term 
was  used  to  represent  interactions  due  to  bond  stretching 

where  Ubo„d  is  the  energy  due  to  bond  stretching,  kb  is  the 
bending  constant  and  rfl  is  the  equilibrium  bond  length. 
Bond  angle  bending  is  also  represented  by  a  harmonic 
potential. 

u^  =  \{8-e,f 

where  Ubend  is  the  energy  due  to  bond  bending,  kg  is  the 
bending  force  constant  and  60  is  the  equilibrium  bond 
angle.  Barriers  to  rotation  about  various  dihedral  angles 
was  controlled  through  a  cosine  series  fit  to  ab  initio 
calculations 

U,or  =XCl(l  +  COS(Z>-£)) 
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where  c,-,  are  torsional  force  constants  and  5  is  the  phase 
angle.  All  force  constants  were  determined  by  fitting 
parameters  to  reproduce  ab  initio  derived  barriers  to  bond 
stretching,  bond  bending  and  rotation  about  dihedral 
angles.  These  calculations  were  performed  in  Gaussian  03 
at  the  HF/6-3  lg+(d,p)  level  of  theory  and  basis  set, 
respectively.  The  ab  initio  calculations  consist  of  scanning 
the  bond  length  or  angle  of  interest,  while  reoptimizing  all 
remaining  degrees  of  freedom.  Equilibrium  bond  lengths 
and  angles  were  determined  through  geometry 
optimization. 


3.  METHODOLOGY 

3.1  Octanol-water  partition  coefficient 

Octanol-water  partition  coefficient  is  the  ratio  of  the 
concentration  of  a  chemical  species  in  octanol  and  in  water 
at  equilibrium.  Octanol  has  been  widely  used  as  a 
surrogate  for  organic  matter  such  as  lipids  since  it  consists 
of  a  nonpolar  tail  group  and  a  polar  head  group.  Octanol- 
water  partition  coefficients  are  of  primary  importance  in 
the  prediction  of  the  bioactivity,  bioaccumulation  and 
transport  of  a  particular  chemical  species.  The 
methodology  used  in  this  work  to  calculate  the  partition 
coefficient  is  a  combination  of  NPT  molecular  dynamics 
and  free  energy  perturbation  (FEP).  The  FEP  method 
involves  slowly  transforming  solute  A  to  solute  B  (either  A 
or  B  is  the  compound  of  interest)  by  scaling  the  interaction 
potential  through 

U(A)  =  AU  +(1-2W 
B  A 


2.2  Bonded  interactions 
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where  A  is  the  scaling  parameter  and  has  values  between 
0  and  1 .  The  Gibbs  free  energy  of  transfer  is  given  by 


AA G  AB  =  AG  B  -  AG  A 
Tr  Tr  Tr 

where  A  Gy,-  A  and  AG/,-  B  are  the  free  energies  for 
transferring  solute  A  and  B  between  the  water  and  the 
water-saturated  octanol  phase  respectively.  Calculation 
of  these  absolute  free  energies  is  difficult.  However, 
because  the  Gibbs  free  energy  is  a  state  function,  the 
relative  Gibbs  free  energies  of  solvation  between  two 
phases  can  be  determined  through  the  construction  of  a 
thennodynamic  cycle,  as  illustrated  in  Fig.  2. 
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Fig.  2.  Thermodynamic  cycle  utilized  to  calculate  octanol- 
water  partition  coefficient. 


The  thermodynamic  cycle  provides  a  means  for 
calculating  the  relative  Gibbs  free  energy  of  transfer 
from  the  following  equation, 

AAG_  AB  =  AG  AB  -  AG  AB 
Tr  oct  aq 

where  AGrpoct)  and  AGrr(w)  are  the  free  energies 
associated  with  the  transformation  of  solute  A  to  solute  B 
in  the  water-saturated  1 -octanol  solution  and  water 
respectively.  The  relative  partition  coefficient  is  now 
expressed  as 

-AA  G  AB 

A  log  K  = - - - 

ow  2.303  RT 

The  absolute  partition  coefficient  of  target  molecule  B  is 
calculated  from  the  reference  molecule  A  by 


log  K  (B)  =  A  log  K  +log  K  ( A ) 

ow  ow  ow 


where  AG  i/wso1  is  the  solvation  free  energy  of  species  i 
in  solvent  water,  pw°  is  the  number  density  of  pure  water 
and  Na  is  the  avagadro's  number.  The  relative  Henry's  law 
constant  is  determined  through  MD  free  energy 
perturbation  from  the  equation, 

Alog10  //  =  AAGml 

where  AAGso1  is  expressed  as 


AA  Gso! 


ags:i(b)-ag:i(a) 

2.303  RT 

AGW(A  — >  B)  -  AGvac(A  — >  B) 
2.303  RT 


using  a  thermodynamic  path  shown  in  Fig.  3. 
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Fig.  3.  Thermodynamic  path  utilized  to  calculate  Henry’s 
law  constant. 


The  absolute  constant  of  B  is  then  calculated  from  A’s 
Henry’s  law  constant  using  equation, 

log  jo  H  ( B )  =  Alog10  H  -  log10  H  (A) 


3.3  Vapor  Liquid  Equilibria  and  Vapor  Pressure 

Gibbs-Duhem  integration  (Kofke,  1993)  is  used  to 
determine  the  phase  coexistence  curve  (temperature  vs 
density)  and  the  vapor  pressure.  With  the  knowledge  of  an 
initial  coexistence  point,  the  Clapeyron  equation  can  be 
integrated  to  provide  an  estimate  of  coexistence  points  at 
other  temperatures.  The  Clapeyron  equation  is  given  by 


3.2  Henry’s  Law  Constant 


Henry's  law  constant  is  the  equilibrium  distribution  of 
a  species  between  gas  and  liquid.  For  dilute  aqueous 
solutions,  it  is  the  ratio  of  the  solute's  partial  pressure  and 
its  aqueous  concentration.  So  it  provides  an  estimate  for 
aqueous  solubility.  Higher  Henry's  law  constant  means 
lower  aqueous  solubility.  Henry's  law  constant  can  be 
expressed  in  terms  of  solvation  free  energy  of  the  solute 
in  water  by  the  following  equation  (Lin  et  al.,  2002) 
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where  P  is  the  pressure,  P  =  I fkT.  Ah  is  the  difference  in 
molar  enthalpies  of  the  coexisting  phases,  Av  is  the 
difference  in  molar  volumes  and  a  indicates  that  the 
derivative  is  taken  along  the  saturation  line.  The  method 
allows  for  the  prediction  of  the  saturation  pressure  at  a 
temperature  AT  away  from  a  known  coexistence  point. 
Given  an  estimate  of  the  saturation  pressure,  NPT 
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Molecular  Dynamics  simulations  are  performed 
simultaneously  for  both  liquid  and  vapor  phase  to 
determine  the  coexistence  densities  and  heat  of 
vaporization.  The  initial  coexistence  point  in  this  work 
was  determined  by  grand  canonical  histogram 
reweighting  Monte  Carlo. 

4.  SIMULATION  DETAILS 


4.1  Octanol-water  Partition  Coefficient 


The  FEP  transformations  employed  in  the  calculations 
are  shown  in  Fig  4. 


Nitrobenzene  MNA 


Fig.  4.  Schematic  of  transformations  used  in  the  FEP 
technique. 

The  reference  solutes  on  the  left  side  of  the  arrow  are 
slowly  transformed  to  the  corresponding  target 
molecules  on  the  right.  The  choice  of  the  reference 
solutes  depends  on  structural  similarity  and  availability 


of  experimental  partition  coefficient  data.  For  each  solute 
of  interest,  two  FEP  simulations  were  performed  at  298  K 
and  1.013  bar;  one  for  the  water  phase,  and  the  other  for 
water-saturated  1-octanol  phase  using  full  periodic 
boundary  conditions.  The  mole  fraction  of  water  in  the 
octanol  phase  was  set  to  the  experimental  value  of  0.255. 
Simulations  were  equilibrated  for  Ins  (nanosecond)  before 
free  energy  calculations  were  initiated.  The  perturbation 
was  carried  out  over  20  windows  where  the  starting  six 
and  the  ending  six  windows  were  unequally  spaced  to 
improve  convergence  at  the  end  points.  This  methodology 
is  known  to  avoid  the  so-called  end-point  catastrophes 
(Beutler  et  al.,  1994;  Pitera  et  al.,  2002)  resulting  from  the 
appearing  and  disappearing  atoms.  After  equilibration, 
FEP  calculations  were  run  for  a  total  of  6  ns  with  100  ps 
(picosecond)  of  equilibration  and  100  ps  of  sampling  for 
each  window  using  a  timestep  of  lfs.  Calculations  were 
performed  in  the  NPT  ensemble  using  the  FEP  module  in 
NAMD.  The  Langevin  piston  Nose-Floover  method 
(Martyna  et  al.,  1994;  Feller  et  al.,  1995)  was  used  to 
control  pressure  and  temperature.  A  non-bonded  cut-off  of 
14  A  is  used  for  all  solutes.  Electrostatic  interactions  were 
calculated  with  the  particle  mesh  Ewald  technique.  The 
TraPPE-UA  (Chen  et  al.,  2001)  force  field  was  used  to 
model  1-octanol  while  the  SPC/E  (Mark  et  al.,  2001)  force 
field  was  used  to  model  water. 

4.2  Henry’s  Law  constant 

The  transformation  and  the  FEP  scheme  used  for 
octanol-water  partition  coefficient  were  also  used  for 
Henry’s  law  constant  calculations.  Two  simulations,  one 
for  the  mutation  in  water  and  one  in  vacuum  were  run  at 
298  K  and  1.013  bar.  Simulations  in  water  used  the  same 
parameters  as  the  FEP  calculations  discussed  in  the 
previous  section.  For  the  calculations  in  vacuum  (in 
vacuo),  an  isolated  hybrid  molecule  was  simulated  without 
periodic  boundary  conditions  and  a  damping  coefficient  of 
10  ps"1  for  Langevin  temperature  control.  The  vacuum  run 
is  carried  out  for  a  total  of  2.4  ns  with  400  ps  of 
equilibration  and  2  ns  of  sampling.  The  in  vacuo 
simulations  require  longer  sampling  times.  A  non-bonded 
cut-off  of  14  A  is  used  for  all  the  solutes. 

4.3  Vapor  Liquid  Equilibria  and  Vapor  Pressure 

One  phase  coexistence  point  near  the  critical  temperature 
was  determined  through  Grand  Canonical  Monte  Carlo 
(GCMC)  simulations  coupled  with  histogram  reweighting 
technique  (Ferrenberg  et  al.,  1988;  Ferrenberg  et  al.,  1989; 
Potoff  et  al.,  1998).  The  insertion  of  molecules  in  the 
GCMC  simulations  was  enhanced  through  multiple  first 
bead  insertions  (Esselink  et  al.,  1995)  and  the  application 
of  the  coupled-decoupled  configurational-bias  monte  carlo 
method  (Martin  et  al.,  1999;  Martin  et  al.,  2006;  Wick  et 
al.,  2000).  The  ratios  of  attempted  moves  were  set  to  60% 
particle  insertions/deletions,  10%  configurational-bias 
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regrowths,  15%  translations  and  15%  rotations.  A  non- 
bonded  cut-off  14  A  was  used.  No  tail  corrections  were 
employed  for  GCMC  calculations.  Treating  this 
coexistence  point  as  the  initial  state,  Gibbs-Duhem 
integration  is  performed  to  generate  the  whole  phase 
diagram.  Given  an  estimate  of  the  saturation  pressure, 
NPT  Molecular  Dynamics  (MD)  simulations  are 
performed  simultaneously  for  both  liquid  and  vapor 
phase  to  determine  the  coexistence  densities,  heat  of 
vaporization  and  subsequently  the  clapeyron  equation  is 
integrated  to  determine  pressure.  The  Langevin  piston 
Nose-Hoover  method  was  used  to  control  pressure  and 
temperature.  A  non-bonded  cut-off  of  14  A  was  used  for 
all  the  simulations  to  maintain  consistency  with  GCMC. 


5.  RESULTS 

5.1  Oetanol-Water  Partition  Coefficients 

The  partial  charges  of  DNAN  and  MNA  determined 
from  Gaussian  03  were  scaled  by  a  factor  of  0.94  to  fit 
the  experimental  octanol-water  partition  coefficients. 
Experimental  octanol-water  partition  coefficients  for  the 
reference  solutes  were  obtained  from  the  literature 
(Schultz,  1999;  Schultz,  1982).  Since  no  experimental 
data  exists  for  either  DNP  or  NTO,  the  same  charge 
scheme  was  applied  to  these  materials.  The  partition 
coefficients  predicted  by  the  FEP  method  are  presented 
in  Table  1  alongside  values  predicted  by  Toghiani  et  al. 
(Toghiani  et  ah,  2008)  through  group  contribution  (Broto 
et  ah,  1984)  and  COSMOthermX  software  and  the 
experimental  data  (Boddu  et  al,  2008;  Boddu  et  al, 
2008).  The  octanol-water  partition  coefficients  predicted 
for  DNAN  and  MNA  are  in  very  good  agreement  with 
the  experimental  values.  The  values  predicted  for  DNP 
and  NTO  compare  well  with  the  group  contribution 
values  but  slightly  deviate  from  COSMOtherm 
predictions. 


TABLE  1.  PREDICTED  OCTANOL-WATER 
PARTITION  COEFFICIENT  (LOG  Kow) 


Compound 

FEP 

Group 

Contribution 

COSMO 

thermX 

Exp 

DNAN 

1.78 

1.70 

1.92 

1.64 

MNA 

2.02 

1.57 

0.80 

2.04 

DNP 

-0.65 

-0.82 

-1.19 

- 

NTO 

1.03 

0.9 

0.37 

- 

5.2  Henry’s  Law  Constants 

Henry’s  law  constants  predicted  for  the  newly 
developed  force  fields  for  DNAN,  MNA,  DNP  and  NTO 
are  shown  in  Table  2.  The  Henry’s  law  constants  for  the 


reference  molecules  were  obtained  from  the  literature 
(Hine  et  ah,  1974).  No  direct  Henry’s  law  constant  has 
been  reported  in  literature  for  pyrazole  although  we  were 
able  to  calculate  that  using  the  relation 

H  =  pi  S 

where  p  is  the  vapor  pressure  and  S  is  the  solubility  in 
water  at  298  K.  The  vapor  pressure  of  pyrazole  at  298  K  is 
3.638  Pa  (Jimenez  et  ah,  1987)  and  the  solubility  in  water 
at  298  K  is  19.4  mohkg  of  water  (Wiley,  1967). 


TABLE  2.  PREDICTED  HENRY’S  LAW  CONSTANTS 
(LOG  H) 


Compound 

FEP 

Exp 

DNAN 

-4.40 

-3.25 

MNA 

-3.84 

-3.60 

DNP 

-6.59 

- 

NTO 

-11.38 

- 

While  DNAN  Henry’s  law  constant  deviates  slightly  from 
the  experiment,  the  value  predicted  for  MNA  agrees  well 
with  the  experiment. 

5.3  Vapor  Liquid  Equilibria  and  Vapor  Pressure 

The  phase  coexistence  curve  and  the  Clausius-Clapeyron 
plot  generated  through  Gibbs-Duhem  integration  are 
shown  in  Fig  5  &  6  respectively. 


Fig.  5.  Vapor  liquid  equilibria  of  DNAN  (circle) 
and  MNA  (square) 

Critical  temperatures  and  densities  are  computed  by  fitting 
the  saturated  liquid  and  vapor  densities  to  the  density 
scaling  law  for  critical  temperature  (Rowlinson  and 
Widom,  1982) 
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Pliq  P vap  =  B(T-Tr) 


and  the  law  of  rectilinear  diameters  (Rowlinson  and 
Swinton,  1982) 


P,-  +P 
liq  vap 


=  p  +  A(T  —  T  ) 
c  c 


where  (3  =  0.325  is  the  critical  exponent  for  Ising-type 
fluids  in  three  dimensions  (Privman  et  al.,  1998)  and  A 
and  B  are  constants  fit  to  simulation  data. 


physical  properties  and  partitioning  of  insensitive 
munitions.  Such  predictive  capability  may  lead  to 
improved  cost  savings  and  reduced  development  time  if 
candidate  molecules  can  be  pre-screened  on  computer 
before  synthesis.  Force  fields  have  been  presented  for 
DNAN,  MNA,  DNP  and  NTO.  The  predictions  of 
simulation  for  the  octanol-water  partition  coefficient  and 
Henry’s  law  constant  were  found  to  be  in  good  agreement 
with  the  spare  experimental  data  available  in  the  open 
literature. 
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Fig.  6.  Clausius-Clapeyron  plot  of  DNAN  (circle) 
and  MNA  (square) 

Normal  boiling  points  and  critical  pressures  were 
determined  from  the  Clausius-Clapeyron  plots  through 
linear-interpolation.  The  critical  parameters  and  the 
normal  boiling  point  are  presented  in  Table  3. 


TABLE  3.  PREDICTED  CRITICAL  PARAMETERS 
AND  BOILING  POINT  FOR  DNAN  AND  MNA 


Compound 

Tc 

(K) 

Pc 

(kg/m3) 

Pc 

(bar) 

Tb 

(K) 

DNAN 

900.01 

40.33 

374.54 

623.75 

MNA 

783.47 

48.45 

345.99 

511.90 

No  experimental  data  exists  for  either  DNAN  or  MNA  to 
make  a  comparison. 

CONCLUSIONS 

In  this  work,  we  have  demonstrated  the  potential  of 
atomistic  computer  simulations  in  the  predictions  of 
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